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Abstract. The Atmospheric Imaging Assembly in the Solar Dynamics Observatory 
provides full Sun images every 12 seconds in each of 7 Extreme Ultraviolet passbands. 
However, for a significant amount of these images, saturation affects their most intense 
core, preventing scientists from a full exploitation of their physical meaning. In 
this paper we describe a mathematical and automatic procedure for the recovery of 
information in the primary saturation region based on a correlation/inversion analysis 
of the diffraction pattern associated to the telescope observations. Further, we suggest 
an interpolation-based method for determining the image background that allows the 
recovery of information also in the region of secondary saturation (blooming). 
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The Solar Dynamics Observatory (SDO) [I] is a solar satellite launched by NASA on 
February 11 2010. The scientific goal of this mission is a better understanding of how 
the solar magnetic field is generated and structured and how solar magnetic energy is 
stored and released into the helio- and geo-sphere, thus influencing space weather. SDO 
contains a suite of three instruments: 

• The Helioseismic and Magnetic Imager (SDO/HMI) [2] has been designed to study 
oscillations and the magnetic field at the solar photosphere. 

• The Atmospheric Imaging Assembly (SDO/AIA) |3] is made of four telescopes, 
providing ten full-Sun images every twelve seconds, twenty four hours a day, seven 
days a week. 

• The Extreme Ultraviolet Variability Experiment (SDO/EVE) [1] measures the 
solar extreme ultraviolet (EUV) irradiance with unprecedented spectral resolution, 
temporal cadence, accuracy, and precision. 

The present paper deals with an important aspect of the image reconstruction 
problem for SDO/AIA jS, 6, 7j. The four telescopes of such instrument capture images 
of the Sun’s atmosphere in ten separate wave bands, seven of which centered at EUV 
wavelengths. Each image is a 4096 x 4096 square array with pixel width in the range 
0.6 — 1.5 arcsec and is acquired according to a standard CCD-based imaging technique. 
In fact, each AIA telescope utilizes a 16-megapixel CCD divided into four 2048 x 2048 
quadrants. As typically happens in this kind of imaging, AIA CCDs are affected by 
primary saturation and blooming, which degrade both quantitatively and qualitatively 
the AIA imaging properties. Primary saturation [E] refers to the condition where a 
set of pixel cells reaches the Full Well Capacity, i.e. these pixels store the maximum 
number possible of photon-induced electrons. At saturation, pixels lose their ability to 
accommodate additional charge, which therefore spreads into neighboring pixels, causing 
either erroneous measurements or second-order saturation. Such spread of charge is 
named blooming [E] and typically shows up as a bright artifact along a privileged 
axis in the image. Figure [l] shows a notable example of combined saturation and 
blooming effects in an SDO/AIA image captured during the September 6, 2011 event. 
The recovery of information in the primary saturation region by means of an inverse 
diffraction procedure is the main goal of the present paper. Further, we also introduce 
here an interpolation approach that allows a robust estimate of the background in the 
diffraction region as well as reasonable estimate of the flux in the central blooming 
region. 

The optical setup of each AIA telescope is characterized by structures with uniform 
wire meshes used to support the thin filters that create the EUV passbands. The 
interaction between the incoming EUV radiation and the grids generates a diffraction 
effect 0 depending on the source intensity Jo- This effect can be easily observed when 
/ 0 is so large that the diffracted flux is comparable to the background level. When 
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Figure 1 . SDO/AIA image of the September 6 2011 event, captured by means of the 
13lA passband at 22:19:09 UT for ~ 2.9 sec exposure duration. The left panel shows 
the position of the explosion on the full disk of the Sun. The zoom in the right panel 
clearly illustrates the presence of primary saturation, blooming and diffraction fringes. 


this situation occurs, it often happens that the incoming flux generates a signal that 
exceeds the saturation level of the AIA CCDs, which is 16383 DN pixel -1 (where DN 
stands for Data Number). This fact has a very interesting mathematical implication: 
all information on the radiation flux which is lost due to primary saturation is actually 
present in the diffraction pattern and therefore the signal in the primary saturation 
region can in principle be restored by solving the inverse diffraction problem (we point 
out that pixels contained in the blooming region generate diffraction effects that are 
negligible with respect to the ones corresponding to the primary saturation region). 

The present paper addresses the de-saturation problem for SDO/AIA having been 
inspired by the heuristic and semi-automatic approaches described in ma and 
respectively. More specifically, we describe here a fully automatic numerical method that 
utilizes correlation and statistical regularization to recover the image information in the 
primary saturation region and applies interpolation to ameliorate the effects of blooming. 
The first pillar of our approach is the definition of a forward model for SDO/AIA 
data formation. In fact the point spread function (PSF) of each passband can be 
approximated as the sum of a core PSF, which is modeled by a two-dimensional Gaussian 
function, and a diffraction PSF, that describes the diffraction pattern corresponding to 
a point source. Therefore the forward model is encoded by the linear integral operator 
whose integral kernel is the sum of the two PSFs. However, the difficult issue here 
is to define the domain of the diffraction PSF, i.e. to automatically segment the 
primary saturation region with respect to the blooming one. We solved this problem 
by combining correlation and thresholding. Once the forward operator is defined, the 
second pillar of our approach performs image reconstruction by means of an Expectation 
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Maximization (EM) algorithm [l_2j applied to the inverse diffraction problem. EM is a 
maximum likelihood technique working when the measured data are affected by Poisson 
noise and the solution to reconstruct is non-negative. We note that AIA data are only 
approximately Poisson, since the system only records the charge and then divides it 
by the average charge per photon to get the Data Number. On the other hand, the 
source distribution in the saturated region, is certainly positive. Further, in this paper 
we will take advantage of a very effective, recently introduced stopping rule for EM, 
which guarantees the right amount of regularization by means of a criterion with solid 
statistical basis im The effectiveness of this de-saturation method is verified against 
synthetic data and by reconstructing the source distribution in the saturated regions of 
SDO/AIA maps of the September 6 2011 event, acquired at different time points. 

The plan of the paper is as follows. Section 2 models the forward problem. Section 
3 describes the image reconstruction method utilized for the solution of the inverse 
diffraction problem. Section 4 performs a numerical validation of the de-saturation 
approach in the case of synthetic data mimicking the SDO/AIA signal formation process. 
An example of how the method works in the case of real data is described in Section 5. 
Finally, our conclusions are offered in Section 6 while an appendix illustrates the way 
we estimate the background and recover information in the blooming region. 

2. The forward problem 

As shown in the Introduction, during many observations an SDO/AIA image presents 
a rather complex structure. Using a lexicographic order for the image pixels and if I 
with size N is one of such images, then in / it is possible to point out five different sets 
of pixels: 

(i) The set of saturated pixels 

S 1 = {i E A/", U = 16383 DN pixel -1 } , (1) 

where N is the set of natural numbers ranging from 1 to N. 

(ii) The subset S C <S" of saturated pixels which are affected by primary saturation. 

(iii) The subset B C <S" of saturated pixels which are affected by blooming. Of course 
we have that 

S' = SUB S O B = 0 . (2) 

(iv) The set of pixels F such that F fl S' = 0 and where the diffraction fringes occur. 

(v) The complement of S' U F (as we will see this set does not play any role in the 
de-saturation process). 

In general, the data formation process in AIA, which generates the image /, is 
the result of the discretization of the convolution between the telescope PSF and the 
incoming photon flux. In a finite-dimension setting, we can introduce the matrix 


A = Ad + Ac , 


( 3 ) 
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which is the sum of the two N x N circulant matrices Ap, associated to the diffraction 
component of the AIA PSF, and Ac, associated to the diffusion component. Therefore 
the image I is given by 

I = Ax = Apx + Acx , (4) 

where x is the vector of dimension N obtained by discretizing the incoming photon 
flux. We now define the sub-matrix Af, : of A D that maps the vector 

x of the values of the photon flux coming just from S onto the vector made of the 
diffraction fringes; here ffS and ffF are the cardinality of S and F, respectively. Since 
the diffraction effects are negligible for pixels outside the core S, we have that the 
diffraction pattern If = {I, , i E F} is approximated by the matrix times vector 
product 

If = A^)X T BGf , (5) 

where BG := Aqx is the total background and BGf '■= ( Aqx)f is its restriction onto 
F. This equation represents the forward model for AIA imaging we are interested in in 
this paper, and is completely defined once S, F and BGf are explicitly estimated. We 
point out that this model neglects the diffraction of the background region on itself as 
well as the diffraction of the bloomed region. It follows that, in this context, the overall 
background BG is the image deprived by the diffraction effects. 

We first observe that F is directly related to S. In fact, if S is known, then 

F = {ieM\S' , (A D xs)i^ 0}, (6) 

where xs is a vector with size N whose components are 1 in A and zero elsewhere. 
Equation (JdJ) points out the set of pixels outside the saturation region S', illuminated 
by the diffraction pattern produced by point sources located in the primary saturation 
region S. On the other hand, if xs 1 is a vector with size N whose components are 1 in 
S' and zero elsewhere, 

F'= {i e AT \ S' , (A DXS ')i^ 0}, (7) 

would point out the set of pixels outside S' illuminated by the diffraction pattern 
produced by point sources located in the overall saturation region S'. Based on the 
observation that the diffraction effects associated to blooming are negligible with respect 
to the ones associated to primary saturation, equations (|6j) and ([7]) suggest that the 
segmentation of S in S' can be obtained by a simple correlation analysis. In fact, let 
us consider first the ideal case where no background affects the fringe data Ip ; then the 
correlation vector is the back-projection 

C = (A s d ) t I f ,, (8) 

where Ip> = {U , i G F'} and {Af ) ) T is the transpose of the matrix Af), i.e. the sub¬ 
matrix of Ad that maps the vector of the photon flux coming from S' on the image 
data in F' (the size of A% is ffF' x if S'). Given C in Q, S would be identified by 
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the positions associated to the components of C larger than 16383 DN pixel 1 . We now 
observe that, if is normalized as 

#F' 

£(4>‘j=l j =!,■■■, (9) 

2=1 

then equation can be interpreted as the first iteration of Expectation-Maximization 
(EM) in absence of background and with initialization given by a unit vector. Therefore, 
in order to account for the presence of background, we generalize the computation of 
the correlation by using 


C (d = c (o) . {A 3 d) t 


__ 

If c^Tbcf 


( 10 ) 


which is the first EM iteration when C'^ 0 -’ is a generic initialization vector, A/ is not 
normalized as in (|9]) and the vector BGpi containing the background values in F' is 
included (we notice that in © and from now on, the symbol • and the fraction should 
be intended as element-wise). It follows that, in order to explicitly compute in (10) 
we need to estimate the background vector BGf' and to select C (0 K In the Appendix 
we show a simple way to estimate the background on the whole image. From it a 
reasonable choice for the initialization is C ((y> = BGs>, he. the background values in S'. 
Once computed C <L> via (10), the primarily saturated region S is given by the positions 
of whose corresponding components are more intense than 16383, the corresponding 
F is given by equation @ and the forward model ([5]) is completely defined. 


3. Inversion method for de-saturation 


The identification of the saturation region by means of correlation allows the definition 
of the SDO/AIA de-saturation problem as the linear inverse problem of determining 
x from the measured data Ip and an estimate BGf of the background in F, when x, 
If and BGp are related by ([5]) . Since the noise affecting the measured data Ip has an 
approximate Poisson nature, the likelihood function, i.e. the probability of obtaining the 
realization Ip of the data random vector given the realization x of the solution random 
vector can be written as 


iZ- Q-( A D x + BG F)i 

p(If\x) = n- TJ-r, - (A s d x + BG F )i ,r) ‘ 


(ii) 


2=1 


A classical statistical approach to the solution of © considers the constrained Maximum 
Likelihood problem 


maxp(Ip\x) 


x > 0. 


( 12 ) 


Expectation Maximization (EM) solves this problem iteratively by means of [T 2j 

If 


X^k+l) __ x 


« . 


S\T 


(A/) 


A s d xW + BG} 


(13) 
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where an appropriate stopping rule introduces a regularization effect. In this application 
we utilize the same statistics-based stopping criterion introduced in ra and successfully 
applied to the reconstruction of RHESSI images [14]. This approach is based on the 
observation that in the case of Poisson noise, standard regularization, referring to a 
specific data vector and concerning the point-wise convergence of a one-parameter family 
of regularizing operators, cannot be applied. Therefore a new definition of asymptotic 
regularization must be introduced, in order to account for the fact that in the Poisson 
case, convergence must hold when the signal-to-noise ratio of the data grows up to 
infinity. More formally: 

Definition 3.1. An operator R : R+ F C R# F —¥ R^ 5 is called asymptotically 
regularizable on a cone C C R+ F if there exists a family of continuous operators 

R {k) : C ->• R* s (14) 

with k G N (or RJ and a parameter choice rule 


k : R_|_ x C —^ Id 


(15) 


such that, for 5 > 0 and Ip, I F £ C 

lim sup{\\RW s ’ i5 f»{I s f )-R(I f )\\ | \\I S F -I F \\<6 , ||/ F || > L } = 0 , (16) 

L—t oo 

where L > 0, I F = 4/II4II? Ip = If/H-TfII, and || • || is the Euclidean norm. Such a 
pair ({R^}, k) is called an asymptotic regularization method for R in C. 


We notice that the EM algorithm can be thought as a family of operators Rf k 1 by 
taking the concatenation of the first k iterations and by thinking it as a function of 
the entry data Ip, i.e. 

R {k \lFi BGp ) = if{i F) BG F ) ° ° i((i f ,bg f ){ 1) (17) 

'-V-" 

k times 


where ip(i F ,BG F ) is th e EM iteration a4 fc+1 ) = f>(i F) BG F ){ x ^) defined in equation (13) 
and 1 the unit vector. Since the EM algorithm is convergent we can define the limit 
operator with no background, i.e 


R(I f ) := lim R {k \l F ,0) (18) 

k—> oo 

for every I F e C [15]. Now we prove that, in the presence of a given background BGp , 
the EM algorithm is an asymptotic regularization for its limit operator R when using 
the following stopping rule [13] . 


Definition 3.2. We call KL-KKT (Kullback Leibler - Karush Khun Tucker) stopping 
rule the function 

k(5, 4) := inf {k e N | P {k \l 6 F , BGp) < rQ^\l s F , BG f )} (19) 


with t > 0 , 


p( k \l F , BGp) := 


x « • ( A s d ) t f 1 


T6 
1 F 


A s d x^ + BGf 


2 


( 20 ) 
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Q*Hl‘ F ,BG F ):=fl<' ( ^ )2(:I< ‘ ))2 

2=1 


A s d x^ + BGi 


( 21 ) 


where division by vectors, (Af>) 2 and (x^) 2 indicate component-wise operations and 
x {k) — x( k \lp, BGf) as in Til 


We now give four technical but easy lemmas. 

Lemma 3.1. Let L > 0 and ^€C. For the EM algorithm fi(i F ,BG F ) defined in equation 
(13) the following relations hold true: 

1p{L-I F ,BG F )(€) = L fj(I F ,BG F )(0 

and 

i>(L-I F ,BG F ){L1i) = L 1p(l F ,BG F /L){€) ■ 

Proof. It follows from computations. 


( 22 ) 

(23) 

□ 


The previous lemma means that, given an initialization, say x^°\ with positive 
components, for k > 2 the k- th EM iteration applied to signal L ■ Ip with background 
BGp is a multiple of the A-tli iterate applied to Ip with background BGp/L, i.e. 
(L • I f ,BGf) = L x^ k \lp, BGp/L). Similar properties hold for the functions P^ 


x 


(k) 


and of Definition 3.2 


Lemma 3.2. For the KL-KKT rule defined above we have 
P {k \L ■ I F , BG f ) = L 2 P < ' k \l F , BGp/L) 

and 

Q (k) (L • I Fl BG f ) = LQ( k \l F , BGp/L) 

for k > 2. 


(24) 

(25) 


Proof It follows by applying equations (24) and (25) to definition (3.2) of P lk) and 

Q™. □ 

The following lemma is a readily generalization of the well-known flux preservation 
condition for the EM algorithm in presence of positive background. 


Lemma 3.3. The EM method x^ k \lp, BGp) as defined in equation 13 has the following 
property 

#F #F 

< YPO, (26) 

2=1 2=1 

and the equality holds only if BGp = 0. 


Proof This follows by replacing k with k + 1 and then by using the explicit form of the 
algorithm (13). □ 
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Lemma 3.4. The functions Q (k> have positive lower and upper bounds for every pair 

C i 5 f,bg f ). 


Proof. By using standard inequalities between norms and Lemma [3~3l we get an upper 
bound for Q^ k \ i.e. 

* F " a ^ 2 (xW(I 5 f ,BG f )) 2 ' 


Q «\r F ,BG F) < ^(%,^ ItiBGr) 


#F #F 

2=1 2=1 


(27) 


and also a lower bound 

Q^ k \l F , BG 


F) > 


T*I 1 ((A s D n^\ii,BG F )f) r 


> 


+ IIbGfIU 

E£(^ b (4,«g f )) 


) 


>o, 


# F Z* F A4)i + \\ bg f \ 

since A^x^ cannot tend to 0 as the KL divergence should tend to infinity. 
Now we can prove the main 


(28) 


□ 


Theorem 3.1. If the EM iteration in equation (13) applied to Ip is stopped at the 
k* iterate with k * := k(S,I F ) according to the KL-KKT criterion (Definition \3.3\ ) then 
fc* < (X) and x^\P F , BG F /\\Ip\\) tends to a solution of problem uM) with data I F (with 


no background) as ||/i?|| —> oo. If, in addition, the solution is unique, then the method 
is an asymptotic regularization for its limit operator (18) on the positive cone M+ F . 


Proof. For the Lemma 3.4 the function Q^ in the definition (19) is bounded away 


from 0. Moreover the function P 00 tends to 0 as k —» oo by construction. Hence, 
k(5, Ip) < oo for each (5, Ip). Let us denote by k* := k(5, Ip) the stopping index. Since, 
by Lemma 3.2, the condition in the definition (19) can be written as 


p( k *\l S p,BGi 


141 


< 


T 


Q (kti a F > bg f /\\i f \\) - ii4ii 


it results that 


lim C'-'di. BG, 


M = o 


(29) 


(30) 


where /) > \lr \ — S, as Q i! '~ 1 is bounded. This means that x^(I F , BG f /\\I f \\) 

satisfies the KKT conditions as \ \I F \\ tends to infinity, which was to be shown. 

If the solution of problem (12) with data (4,0) is unique, say x, then x = R(I F ) 
where R(I F ) is the limit operator defined in (18), and x^ k *\lp, BG F /\\Ip\\) —> x as 
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peak value [10 4 ] 

a [arcsec] 

position [arcsec] 

4 

1.5 

(-20,-5) 

3 

2.5 

(5,5) 

5 

1.0 

(-2,7) 


Table 1 . The physical and geometric parameters associated to the synthetic sources 
used in the first validation test for the inverse diffraction approach. 


We point out that EM in (13), together with the KL-KKT stopping rule given in 
Definition |3.2[ provides a regularized vector x in the object space of the reconstructed 
images, while our aim is data de-saturation in the data space, where saturation 
actually occurs. Therefore our procedure ends with the projection of the KL-KKT 
EM reconstructed x in order to construct the desaturated image Idesat- The scheme is: 


Idesat < 


1 C 


A s r 

in S 


BG b 

in B 

(31) 

Ip — AfjX 

in F 

I 

V 

in JV \ ( S' U F ) . 



is the sub-matrix of Ac 

that maps the vector of the 


values of a; in S' onto the corresponding vector in S and BGp is the restriction to the 
blooming region B of the background estimated as described in the Appendix. 


4. Numerical validation 

In order to validate the effectiveness of this approach for restoring information in 
the primary saturation region we consider two simulations that mimic the presence 
of primary saturation in SDO/AIA data with two different levels of adherence with 
experimental conditions (the way blooming can be eliminated from experimental images 
is illustrated in the Appendix and some examples are given in the next Section). 

In the first simulation, the configuration that generates the synthetic data is 
represented by three two-dimensional symmetric Gaussian functions whose geometrical 
characteristics are described in Table Q] We added a constant offset to this simulated 
configuration in order to mimic the presence of background and numerically convolved 
the resulting image with the global AIA PSF A in equation (|3]). The resulting blurred 
and diffracted image (see Figure [2j top left panel) was affected by Poisson noise and the 
result was artificially saturated by setting to 16383 DN pixel -1 all grey levels greater than 
this value (see Figure [2j top right panel). We finally applied our numerical automatic 
procedure described in Section 3 and obtained the result in Figure [2| bottom left panel. 
Finally, the bottom right panel shows the level of accuracy with which the reconstruction 
method is able to recover information inside the saturated region of the image (the root 
mean square error in the saturation region is of the order of 9%). 
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Figure 2. Validation of the correlation/inversion method against synthetic data. 
Top left panel: the ground-truth image obtained from the synthetic configuration by 
applying diffraction and blurring. Top right panel: artificially saturated image with 
Poisson noise added. Bottom left panel: de-saturated image. Bottom right panel: 
lexicographic representation of the grey-level values for the pixels in the saturated 
region: ground truth values (dashed line); saturated values (solid line); reconstructed 
values (black line). 


The second simulation created a synthetic saturated dataset starting from the non- 
satnrated real AIA 4096 x 4096 image in Figure |3j top left panel. The simulation process 
is implemented by means of the following steps: 

(i) The diffraction fringes and the diffusion blurring are eliminated by applying a 
deconvolution step based on EM (with the global PSF), which provides the map in 
Figure |3j top right panel. 

(ii) The pixels in the brighter part of the image have been carried over the saturation 
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threshold, to obtain the image in Figure |3j middle left panel. To do this we applied 
the following re-scaling procedure: denoting with x the pixel intensity before the 
rescaling, the rescaled pixel intensity is 


%rescaled 


mM-x* ~ r M(l-m) ~ \ 

M—x* ^ ' M-x* — 

X X < X* 


(32) 


where x* = 0.25 M, M is the maximum intensity in the image and m — 12 in the 
figure. 

(iii) The image in Figure [3j middle left panel, is first convolved with the core PSF Ac 
to construct the prototype, in Figure [3j middle right panel, of the ideal image that 
would be recorded by AIA if no diffraction and no saturation occurred. This is the 
ground-truth we want to restore. 

(iv) The same image in Figure [3| middle left panel, is now convolved with the global 
PSF A and the result is first affected by Poisson noise and then saturated by setting 
to 16383 DN pixel -1 all the pixel values over the saturation threshold. The result 
of this step is in Figure [3j bottom left panel. 

(v) Finally, we apply the de-saturation method described in Section 3 to this image, to 
obtain the reconstruction in Figure |3| bottom right panel. 


In order to quantitatively assess the reliability of this procedure we have computed, for 
different values of m, the C-statistics 


Cstatib, x) = — y2(I F )i log 


(l 


F i 


+(ApX+BGF)i—(lF)i, (33) 


(A s d x + BGp)i 

that measures, according to the Kullbach-Leibler topology, the discrepancy between the 
data in the region of the diffraction fringes and the expectation corresponding to the 
reconstructed source. Table [2] shows such C sta t values; computes the root mean square 
errors in S between the ground truth and the de-saturated images; and, finally, compares 
the sum Tp of the values of the original flux Ip from the pixels in the diffraction fringes, 
with the sum T|7 ed of the values of the flux in the diffraction fringes predicted by the 
reconstructed source, again for different values of m in equation (32). 


5. Application to real data 

The processing and, specifically, the de-saturation of SDO/AIA data is truly a big 
data issue: indeed, AIA data include images of the Sun in 7 EUV wavelengths every 
12 seconds since February 2010. Here we just focused on a small set of examples 
and, in particular, we considered the single event occurred on September 6 2011, 
detected at the four different wavelengths 94A, 13L4, 171A, and 193A, around the 
same acquisition time, i.e. 22:18:50 UT, 22:19:25 UT, 22:16:48 UT, and 22:16:43 UT, 
respectively. Figure [4] compares the original saturated images with the reconstructed 
ones. We used the correlation/inversion process described in Section 3 in order to 
recover the information lost due to primary saturation and the interpolation procedure 
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Figure 3. The de-saturation method at work in the case of an experimental image 
synthetically saturated. Top left panel: the original image with highlighted the region 
to de-convolve. Top right panel: the deconvolved image. Middle left panel: re-scaled 
image. Middle right panel: prototype of the ideal un-diffracted, un-saturated image. 
Bottom left panel: saturated image affected by Poisson noise. Bottom right panel: de- 
saturated image. The images to compare are this last one and the one in the middle 
right panel. 
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m 

C'stat 

r F (io 6 ) 

T|7 ed (10 6 ) 

RMS (%) 

6 

4.66 

3.18 

3.26 

4.36 

9 

3.11 

4.33 

4.43 

1.72 

12 

2.29 

5.99 

6.06 

0.95 

15 

2.08 

6.78 

6.82 

0.56 

18 

1.67 

7.51 

7.52 

0.36 


Table 2. The de-saturation method at work in the case of an experimental 
image synthetically saturated: C-statistic, original and reconstructed flux intensity 
in correspondence with the diffraction fringes, and root mean square error for several 
re-scaling intensities. 


wavelength (A) 

time (UT) 

7>(10 6 ) 

T| red (10 6 ) 

Cstat 

94 

22:18:50 

3.44 

3.50 

2.83 

131 

22:19:25 

4.51 

4.55 

1.97 

171 

22:16:48 

6.73 

6.76 

3.66 

193 

22:16:43 

1.16 

1.15 

8.97 


Table 3. The de-saturation method at work in the case of an experimental 
image recording during the September 6 2011 event. The fluxes are computed in 
correspondence with the diffraction fringes. 


described in the Appendix to reduce the effect of blooming. The figure visually 
demonstrates the effectiveness of the de-saturation process. However, in order to provide 
a quantitative assessment of such effectiveness, in Table [3] we provide the C-statistic 
values in correspondence of the diffraction fringes for the four different wavelengths. 

Finally, Figure [5] shows the effectiveness of the method in the case of data recorded 
at a time point when the saturation effect was really dramatic. However the de¬ 
saturation method shows once more its power in both the reconstruction of the primary 
information and in the data interpolation for the blooming region. It is interesting to 
note that, for this case, the C-statistic value is higher (~ 14.03) than for the cases 
described in Table [3} This is most likely due to the fact that our algorithm does not 
account for the wavelength-dependent dispersion of the AIA PSFs (the formulation of 
a multi-wavelength approach to this de-saturation problem is in progress). 

6. Conclusions 

SDO/AIA images are strongly affected by both primary saturation and blooming, 
that may occur at all different wavelengths and acquisition times, even in the case 
of flaring events characterized by a moderate peak flux. This paper describes the first 
mathematical description of a robust method for the de-saturation of such images at 
both a primary and secondary (blooming) level. The method relies on the description 
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Figure 4. The de-saturation method at work for experimental data recorded during 
the September 6, 2011 event. First row panels: experimental and de-saturated 
images for the 94 A bandwidth at 22:18:50 UT. Second row panels: experimental 
and de-saturated images for the 131A bandwidth at 22:19:25 UT. Third row panels: 
experimental and de-saturated images for the 111 A bandwidth at 22:16:48 UT. 
Fourth row panels: experimental and de-saturated images for the 193 A bandwidth 
at 22:16:43 UT. 
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Figure 5. The case of the data collected at 22:19:09 UT, on September 6 2011, by 
means of the 13lA passband. In this case the saturation effects are really impressive 
but the method is still able to recover the information in the saturated region. 


of de-saturation in terms of inverse diffraction and utilizes correlation and Expectation- 
Maximization for the recovery of information in the primarily saturated region. This 
approach requires to compute a reliable estimate of the image background which, for 
this paper, has been obtained by means of interpolation in the Fourier space. The 
knowledge of the background permits to recover information in the blooming region in 
a very natural way. 

The availability of an automatic procedure for image de-saturation in the SDO/AIA 
framework may potentially change the extent with which EUV information from the Sun 
can be exploited. In fact, armed with our computational approach, many novel problems 
can be addressed in SDO/AIA imaging. For example, one can study the impact of the 
choice of the model for the diffraction PSF on the quality of the de-saturation. In this 
paper we used a synthetic estimate of the diffraction PSF provided by Solar Software 
(SSW) but other empirical or semi-empirical forms can be adopted. Furthermore, this 
technique can be extended to account for the dependance of the PSF from the passband 
wavelengths. Finally, the routine implementing this approach is fully automated and 
this allows the systematic analysis of many events recorded by AIA and their integration 
with data provided by other missions such as RHESSI [T6]J or, in the near future STIX 

nzt 
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The SDO/AIA hardware is equipped with a feedback system that reacts to saturation 
in correspondence of intense emissions by reducing the exposure time. As a result, for 
a typical AIA acquisition along a time range of some minutes, during which saturation 
occurs, the telescope always provides some unsaturated frames that can be utilized to 
estimate the background. A possible way to realize such an estimate is based on the 
following scheme. Let us denote with A and / 3 two unsaturated images acquired at times 
A and A, respectively and with I> a saturated image acquired at A with A < A < 0 
(note that A, A, and / 3 are normalized at the same exposure time). The algorithm for 
the estimate of the background is: 

(i) A and / 3 are deconvolved with EM (using the global PSF) to obtain the 
reconstructions x\ and 7 3 (the KL-KKT rule can be used to stop the iterations). 

(ii) Both X\ and x : > are Fourier transformed by means of a standard FFT-based 
procedure to obtain x\ and x 3 . 

(iii) A low-pass filter is applied to both x\ and x 3 to obtain x 1 and x 3 , respectively. 

(iv) For each corresponding pair of pixels in x{ and x 3 that are not negligible, an 

interpolation routine is applied, both for the real and imaginary part. This provides 
x 2 in correspondence of t — t 2 - 

(v) The resulting vector x nt is Fourier inverted to obtain the interpolated 
reconstruction x™ 1 . 

(vi) The core PSF Ac is finally applied to x™ 1 to obtain I 2 nt in the image domain. 

I 2 nt is a reliable estimate, for time t 2 ? of the background BG = Aqx introduced in 

Section 2. On the other hand, a reliable estimate of the image in the bloomed region is 
provided by the restriction BGb of I 2 nt onto B determined as in Section 2. We finally 
observe that, in this algorithm, the interpolation step is applied in the Fourier domain 
because, after filtering, a lot of pixels are negligible and therefore the computational 
burden of the procedure is notably decreased. 
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